#include "Computations.h"

#include <stdlib.h> 

#include "../Utils/Utils.h"
#include "../SparceMatrix/sparse_matrix.h"
#include "../Containers/NumList.h"
#include "../LinearComputation/Vector.h"
#include "../LinearComputation/Matrix.h"
#include "../LinearComputation/ModularityMatrix.h"



elem* ModularityMult(const sparse_matrix_lst* adjMatrix, const elem* degVector,
		elem MInvers, const elem* fVector, elem norm, const elem* vector)
{
	elem* ret = NULL;
	elem* temp = NULL;
	unsigned size = adjMatrix->n;
	
	NAssert(MemoryAssert(2, 
		ret = (elem*)malloc(sizeof(elem) * size),
		temp = (elem*)malloc(sizeof(elem) * size)));
	
	mult_sparse_list(adjMatrix, vector, ret);
	
	MultVectorByScalar(degVector, ScalarProduct(degVector, vector, size) * MInvers, temp, size);
	SumVectors(ret, temp, ret, size, FALSE);
	
	if (fVector != NULL)
	{
		MultedVector(fVector, vector, temp, size);
		SumVectors(ret, temp, ret, size, FALSE);
	}
	
	if (norm)
	{
		MultVectorByScalar(vector, norm, temp, size);
		SumVectors(ret, temp, ret, size, TRUE);
	}
	
	free(temp);
	
	return ret;
}

EigenpairAndMatrix* PowerMethod(const sparse_matrix_lst* adjMatrix,
		const elem* degVector, const elem* modularityMat, const elem* fVec, 
		elem precision, elem M, unsigned size)
{
	int i = 0;
	EigenpairAndMatrix* ret = NULL;;
	EigenpairAndMatrix* temp = NULL;
	elem norm = 0;
	elem* iterVector = NULL;

	NAssert(MemoryAssert(2,
			ret = (EigenpairAndMatrix*)malloc(sizeof(EigenpairAndMatrix)),
			iterVector = (elem*)malloc(sizeof(elem) * size)));

	norm = Norm(modularityMat, size);
	
	for (i = 1; i < size; i++)
		iterVector[i] = 0;
	iterVector[0] = 1;	
	
	if (NULL != (temp = PowerMethodHelper(adjMatrix, degVector,
			1.0 / M, fVec, norm, iterVector, precision * precision)))
	{
		ret->Vector = temp->Vector;
		ret->Value = temp->Value - norm;
	}
	else
	{
		free(ret);
		ret = NULL;
	}
	
	free(temp);

	return ret;
}


EigenpairAndMatrix* PowerMethodHelper(const sparse_matrix_lst* adjMatrix, 
		const elem* degVector, elem mInvers, const elem* fVector, 
		elem norm, elem* iterVector, elem precision)
{
	EigenpairAndMatrix* ret;
	elem* nextIterVector = NULL;
	elem* temp = NULL;
	elem eigenvalue = 0;
	unsigned size = adjMatrix->n;
	
	NAssert(MemoryAssert(2, 
			nextIterVector = ModularityMult(
					adjMatrix, degVector, mInvers, fVector, norm, iterVector),
			temp = (elem*)malloc(sizeof(elem) * size)));
	
	eigenvalue = ScalarProduct(iterVector, nextIterVector, size) /
				ScalarProduct(iterVector, iterVector, size);
	
	MultVectorByScalar(iterVector, eigenvalue, iterVector, size);
	SumVectors(nextIterVector, iterVector, temp, size, FALSE);
	
	if (ScalarProduct(temp, temp, size) < precision)
	{
		MemoryAssert(3, NULL, iterVector, temp);
		NAssert(MemoryAssert(1,
				ret =(EigenpairAndMatrix*)malloc(sizeof(EigenpairAndMatrix))));
		
		ret->Value = eigenvalue;
		ret->Vector = nextIterVector;
		return ret;
	}
	else
	{
		MemoryAssert(3, NULL, iterVector, temp);
		
		MultVectorByScalar(nextIterVector, 1 / NormVector(nextIterVector, size),
				nextIterVector, size);
		return PowerMethodHelper(adjMatrix, degVector, mInvers, fVector, 
				norm, nextIterVector, precision);
	}
}

elem* SpectricalDiv (const sparse_matrix_lst* adjMatrix, 
		const elem* degVec, const elem* modularityMat, const elem* fVec, 
		elem precision, elem M, unsigned size, elem* globalDeltaQ)
{
	EigenpairAndMatrix* powerProduct = NULL;
	elem* tempVec= NULL, *ret = NULL, deltaQ;
	unsigned i = 0;
	
	NAssert(MemoryAssert(1, 
			powerProduct = PowerMethod(adjMatrix, degVec, 
					modularityMat, fVec, precision, M, size)));
	
	BrutalRoundVector(powerProduct->Vector, size);
	
	NAssert(MemoryAssert(3, 
			tempVec = MultVectorByMatrix(modularityMat, powerProduct->Vector, size), 
			powerProduct->Vector, powerProduct));
	
	deltaQ = 0.5 * ScalarProduct(tempVec, powerProduct->Vector, size);
	
	if(powerProduct->Value > 0 && deltaQ > 0)
	{
		*globalDeltaQ += deltaQ;
		ret = powerProduct->Vector;
	}
	else
	{
		for (i = 0; i < size; i++)
			ret[i] = 1;
		free(powerProduct->Vector);
	}

	MemoryAssert(3, NULL, powerProduct, tempVec);
	
	return ret;
}
